library(ggplot2); library(plotly); library(data.table); library(sf)
## Warning in register(): Can't find generic `scale_type` in package ggplot2 to
## register S3 method.
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
Innovasea uses an Azimuthal Equidistant (AEQD) projection, so going to use that directly rather than switching back to lat/long (WGS84). Filtering to keep only the 90% of detctions with the lowest horizontal positioning error (HPE).
vps <- fread('data/raw/embargo/2021 vps/results/animal/all.csv', col.names = tolower)
vps <- st_as_sf(vps,
coords = c('x', 'y', 'z'),
crs = '+proj=aeqd +lat_0=38.542096 +lon_0=-75.761724 +x_0=1800 +y_0=1800',
remove = F)
vps <- setDT(vps)[hpe < quantile(hpe, 0.9)]
head(vps)
Bring in the lengths at tagging. Most detections are the reference tag, 65011, so the number of actual animal positions we have it much lower than it seems.
tag_data <- fread('data/raw/embargo/sturgeon capture data.csv', col.names = tolower)
vps <- vps[tag_data, on = 'fullid == transmitternumber', nomatch = 0]
head(vps)
nan <- st_read('data/raw/geo/NHD_H_0208_HU4_GDB.gdb',
query = "SELECT OGR_GEOM_WKT AS wkt
FROM wbdhu10
WHERE States LIKE 'D%'")
## Reading query `SELECT OGR_GEOM_WKT AS wkt
## FROM wbdhu10
## WHERE States LIKE 'D%'' from data source `C:\Users\darpa2\Analysis\MarshyhopeAS-2019-22\Data\Raw\GEO\NHD_H_0208_HU4_GDB.gdb'
## using driver `OpenFileGDB'
## Simple feature collection with 7 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -75.97489 ymin: 38.12779 xmax: -75.22449 ymax: 38.99471
## Geodetic CRS: NAD83
nan <- st_read('data/raw/geo/NHD_H_0208_HU4_GDB.gdb',
layer = 'nhdarea',
wkt_filter = nan$wkt)
## Reading layer `nhdarea' from data source
## `C:\Users\darpa2\Analysis\MarshyhopeAS-2019-22\Data\Raw\GEO\NHD_H_0208_HU4_GDB.gdb'
## using driver `OpenFileGDB'
## Simple feature collection with 1 feature and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XYZ
## Bounding box: xmin: -75.97139 ymin: 38.2467 xmax: -75.55369 ymax: 38.70698
## z_range: zmin: 0 zmax: 0
## Geodetic CRS: NAD83
nan <- st_transform(nan, '+proj=aeqd +lat_0=38.542096 +lon_0=-75.761724 +x_0=1800 +y_0=1800')
seg <- st_crop(nan,
xmin = 1300, ymin = 1000,
xmax = 2350, ymax = 2450)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
ggplot() +
geom_sf(data = seg) +
geom_sf(data = vps, aes(geometry = geometry), alpha = 0.3)
## F-ARIS deployment period, by hour
aris <- data.table(
start = seq.POSIXt(
as.POSIXct('2021-08-19 14:00', tz = 'America/New_York'),
as.POSIXct('2021-10-18 11:00', tz = 'America/New_York'),
by = 'hour'
)
)
## Make a dummy "end" column
aris[, end := start]
## Periods where the camera failed
camera_failure <- data.table(
start = as.POSIXct(c('2021-08-22 12:00', '2021-08-26 10:00', '2021-09-07 18:00',
'2021-10-07 04:00', '2021-10-11 11:00'), tz = 'America/New_York'),
end = as.POSIXct(c('2021-08-23 16:00', '2021-09-07 14:00', '2021-09-16 13:00',
'2021-10-07 12:00', '2021-10-14 12:00'), tz = 'America/New_York')
)
## Find which hours of deployment had camera failure
setkey(camera_failure, start, end)
aris <- foverlaps(aris, camera_failure)
## Camera was on if the hour was not during a period of camera failure (i.e., ==NA)
aris[, camera_status := fifelse(is.na(start), 'Camera on', 'Camera off')]
## Remove unneeded columns and rename column that lists the hours
aris[, ':='(start = NULL, end = NULL, i.start = NULL, i.end = NULL,
'hrs' = i.end)]
aris_loc <- st_sfc(st_point(c(-75.76290, 38.54498)), crs = 4326)
aris_loc <- st_transform(aris_loc, '+proj=aeqd +lat_0=38.542096 +lon_0=-75.761724 +x_0=1800 +y_0=1800')
vps_aris <- vps[, dist := as.numeric(st_distance(geometry, aris_loc))]
vps_aris <- vps_aris[dist <= 100]
vps_aris[, .N, by = 'fullid']
ggplot() +
geom_sf(data = seg) +
geom_sf(data = vps_aris, aes(color = fullid, geometry = geometry), alpha = 0.3,
show.legend = F) +
geom_sf(data = aris_loc, shape = 10, size = 5, color = 'red') +
scale_color_viridis_d()
setattr(aris$hrs, 'tzone', 'UTC')
aris[, end := shift(hrs, -1)]
aris <- aris[!is.na(end)]
setkey(aris, hrs, end)
vps_aris[, dummy.end := time]
vps_aris[, time.local := time]
setattr(vps_aris$time.local, 'tzone', 'America/New_York')
vps_aris <- foverlaps(vps_aris, aris, by.x = c('time', 'dummy.end'))[camera_status == 'Camera on']
vps_aris[, .N, by = 'fullid']
ggplot() +
geom_sf(data = seg) +
geom_sf(data = vps_aris, aes(geometry = geometry, color = fullid),
show.legend = F) +
geom_sf(data = aris_loc, shape = 10, size = 5) +
scale_color_viridis_d()
ggplot() +
geom_sf(data = seg) +
geom_sf(data = vps_aris, aes(geometry = geometry, color = fullid),
show.legend = F) +
geom_sf(data = aris_loc, shape = 10, size = 5, color = 'red') +
coord_sf(xlim = c(1590, 1750), ylim = c(2000, 2200)) +
scale_color_viridis_d()
vps_aris[, .(fullid, totallength, time, time.local, datecaptured, dist)]
ggplotly(
ggplot() +
geom_sf(data = seg) +
geom_sf(data = vps_aris, aes(geometry = geometry, color = fullid,
label = totallength,
label2 = time,
label3 = dist),
size = 1,
show.legend = F) +
geom_sf(data = aris_loc, shape = 10, size = 5, color = 'red') +
coord_sf(xlim = c(1590, 1750), ylim = c(2000, 2200)) +
scale_color_viridis_d()
)
## Warning: Ignoring unknown aesthetics: label, label2, label3